#include "mt3d_cuda.h"

void mt3d::save_model(std::string filename, bool debug)
{
    gctl::gmshio modelio;

    std::clog << "Writing Gmsh file ... ";
    modelio.init_file(filename, gctl::Output);
    modelio.set_packed(gctl::NotPacked, gctl::Output);
    modelio.save_mesh(mesh_elems_, mesh_nodes_);
    modelio.save_data("Node tag", node_markers_, gctl::NodeData);
    modelio.save_data("Element tag", elem_markers_, gctl::ElemData);

    if (debug) // save electronic field
    {
        gctl::array<double> node_vol(node_num_, 0.0);
        gctl::array<double> te_xval_real(node_num_, 0.0), te_xval_imag(node_num_, 0.0);
        gctl::array<double> te_yval_real(node_num_, 0.0), te_yval_imag(node_num_, 0.0);
        gctl::array<double> te_zval_real(node_num_, 0.0), te_zval_imag(node_num_, 0.0);
        gctl::array<double> tm_xval_real(node_num_, 0.0), tm_xval_imag(node_num_, 0.0);
        gctl::array<double> tm_yval_real(node_num_, 0.0), tm_yval_imag(node_num_, 0.0);
        gctl::array<double> tm_zval_real(node_num_, 0.0), tm_zval_imag(node_num_, 0.0);

        gctl::vertex3dc* tmp_vert[4];
        double x, y, z, ux, uy, uz, dummy;
        cuDoubleComplex te_fx, te_fy, te_fz, tm_fx, tm_fy, tm_fz;
        efem_sf::ordertype_e o_type;

        double tmp_vol, elem_cndt;
        for (size_t e = 0; e < elem_num_; e++)
        {
            // 双重加权平均 同时考虑三角形面积与电导率的权重 减少空气层界面的畸变
            get_element_conductivity(elem_markers_[e], elem_cndt);
            tmp_vol = mesh_elems_[e].volume() * elem_cndt;

            for (size_t n = 0; n < 4; n++)
            {
                node_vol[mesh_elems_[e].vert[n]->id] += tmp_vol;
            }
        }

        for (size_t e = 0; e < elem_num_; e++)
        {
            tmp_vert[0] = mesh_elems_[e].vert[0];
            tmp_vert[1] = mesh_elems_[e].vert[1];
            tmp_vert[2] = mesh_elems_[e].vert[2];
            tmp_vert[3] = mesh_elems_[e].vert[3];

            get_element_conductivity(elem_markers_[e], elem_cndt);
            tmp_vol = mesh_elems_[e].volume() * elem_cndt;

            for (size_t n = 0; n < 4; n++)
            {
                x = mesh_elems_[e].vert[n]->x;
                y = mesh_elems_[e].vert[n]->y;
                z = mesh_elems_[e].vert[n]->z;

                te_fx = te_fy = te_fz = tm_fx = tm_fy = tm_fz = _zero;
                for (size_t i = 0; i < 6; i++)
                {
                    if (elem_edge_reversed_[e][i]) o_type = efem_sf::Reverse;
                    else o_type = efem_sf::Normal;

                    elsf_.tetrahedron(x, y, z, 
                        tmp_vert[0]->x, tmp_vert[1]->x, tmp_vert[2]->x, tmp_vert[3]->x, 
                        tmp_vert[0]->y, tmp_vert[1]->y, tmp_vert[2]->y, tmp_vert[3]->y, 
                        tmp_vert[0]->z, tmp_vert[1]->z, tmp_vert[2]->z, tmp_vert[3]->z, 
                        i, efem_sf::Value, o_type, ux, uy, uz, dummy, dummy, dummy);

                    te_fx = clcg_Zsum(te_fx, clcg_Zscale(ux, h_X_TE_[elem_edge_index_[e][i]]));
                    te_fy = clcg_Zsum(te_fy, clcg_Zscale(uy, h_X_TE_[elem_edge_index_[e][i]]));
                    te_fz = clcg_Zsum(te_fz, clcg_Zscale(uz, h_X_TE_[elem_edge_index_[e][i]]));
                    tm_fx = clcg_Zsum(tm_fx, clcg_Zscale(ux, h_X_TM_[elem_edge_index_[e][i]]));
                    tm_fy = clcg_Zsum(tm_fy, clcg_Zscale(uy, h_X_TM_[elem_edge_index_[e][i]]));
                    tm_fz = clcg_Zsum(tm_fz, clcg_Zscale(uz, h_X_TM_[elem_edge_index_[e][i]]));
                }

                te_xval_real[mesh_elems_[e].vert[n]->id] += te_fx.x*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_xval_imag[mesh_elems_[e].vert[n]->id] += te_fx.y*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_yval_real[mesh_elems_[e].vert[n]->id] += te_fy.x*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_yval_imag[mesh_elems_[e].vert[n]->id] += te_fy.y*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_zval_real[mesh_elems_[e].vert[n]->id] += te_fz.x*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                te_zval_imag[mesh_elems_[e].vert[n]->id] += te_fz.y*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];

                tm_xval_real[mesh_elems_[e].vert[n]->id] += tm_fx.x*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_xval_imag[mesh_elems_[e].vert[n]->id] += tm_fx.y*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_yval_real[mesh_elems_[e].vert[n]->id] += tm_fy.x*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_yval_imag[mesh_elems_[e].vert[n]->id] += tm_fy.y*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_zval_real[mesh_elems_[e].vert[n]->id] += tm_fz.x*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
                tm_zval_imag[mesh_elems_[e].vert[n]->id] += tm_fz.y*tmp_vol/node_vol[mesh_elems_[e].vert[n]->id];
            }
        }

        modelio.save_data("TE_Ey_real", te_yval_real, gctl::NodeData);
        modelio.save_data("TE_Ey_imag", te_yval_imag, gctl::NodeData);
        modelio.save_data("TM_Ex_real", tm_xval_real, gctl::NodeData);
        modelio.save_data("TM_Ex_imag", tm_xval_imag, gctl::NodeData);
        
    }

    std::clog << "done\n";
    return;
}
